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with Non-homogeneous Boundary Conditions and its Numerical 
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Abstract 

To overcome the weakness of a total variation based model for image restoration, various high order 
(typically second order) regularization models have been proposed and studied recently. In this paper we 
analyze and test a fractional-order derivative based total a-order variation model, which can outperform 
the currently popular high order regularization models. There exist several previous works using total 
a-order variations for image restoration; however first no analysis is done yet and second all tested 
formulations, differing from each other, utilize the zero Dirichlet boundary conditions which are not 
realistic (while non-zero boundary conditions violate definitions of fractional-order derivatives). 

This paper first reviews some results of fractional-order derivatives and then analyzes the theoretical 
properties of the proposed total a-order variational model rigourously. It then develops four algorithms 
for solving the variational problem, one based on the variational Split-Bregman idea and three based 
on direct solution of the discretise-optimization problem. Numerical experiments show that, in terms of 
restoration quality and solution efficiency, the proposed model can produce highly competitive results, 
for smooth images, to two established high order models: the mean curvature and the total generalized 
variation. 

Keywords. Fractional-order derivatives; Total a-order variation; PDF; Image Denoising; Image 
inverse problems; Optimization methods. AMS. 62H35, 65N22, 65N55, 74G65, 74G75 


1 Introduction 

This paper presents a fractional-order derivative based regularizer for variational image restoration. It may be 
used for other imaging models such as image registration. Denote an observed image hy z = z{x), x G ft 
where ft is the bounded domain of the image with d space dimension and has a Lipschitz boundary. Here we 
consider d = 2 and mainly the image denoising problem with an additive noise i.e. assume z = u + rjQ with 
? 7 o representing some unknown Gaussian noise of mean zero and deviation cr, but most results are applicable 
to d > 2 and other noise models. 

1.1 Image inverse problem 

Restoring the unknown u (without any restrictions) from z is an inverse problem. According to the maximum 
likelihood principle [39], most image processing problems involve solving the least-square problem 

min / \P{u) — z\‘^dx, (1) 

“ Ja 

measuring the fidelity to z. For example, P{u) = u for image denoising, P{u) takes the template image 
T(x u(x)) (and z = P(x) for a reference image) for image registration, and P(u) = Pn^(u(x)) for image 
inpainting with fti C ft the subdomain with missing data. 
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The problem 0 is in general ill-posed due to non-uniqueness, therefore how to effectively solve it becomes 
a fundamental task in image sciences. The most popular idea is to regularize it so that the resulting well- 
posed problem admits an unique solution. The classical regularization technique by Tikhonov et. al izg is to 
add a smoothing regularization term into the energy functional to derive the following minimization problem 

min [ \P{u) — z\‘^dx + A / \\/u\'^dx, (2) 

“ Jn Jn 

where A is a positive constant. This model cannot preserve image edges, though it is simple to use. The 
total variation (TV) model by Rudin-Osher-Fatemi [57j or the ROF model 

min / iVnldx, / \P{u) — z\‘^dx = , P{u) = u (3) 

“ Jn Jn 

is widely used, where a is an estimate of the error rjQ between the noisy image z and the true data u. The 
ROF model preserves the image edges by seeking solutions of piecewise constant functions in the space of 
bounded variation functions (BV). A variety of methods based on the TV regularization have been developed 
to deal with the imaging problems such as image restoration [n El uni i 82 ], image registration gg EH | 62 ], 
image decomposition ISI1IM1I31], image inpainting and image segmentation [HITT]. Restoring 

smooth images in some applications where edges are not the main features presents difficulties for the ROF 
model as it can yield the so-called blocky (staircase) effects. Another disadvantage of the model is to the loss 
of image contrasts [52j . It should be remarked that the recently popular method by the iterative regularization 
technique |60] can reduce the staircasing effect and improve on the image contrast to some extent; besides it 
provides a fast implementation. 

1.2 High-order regularization 

To remedy the above mentioned two drawbacks (stairicasing and contrast), two types of alternative regularizer 
to the TV have been proposed in the literature. The first type introduces higher order regularization into 
image variational models EauaiiiisiiiiiEiiiigiHi]. The mean curvature-based variation denoising model 
was studied in [sgEginiiHi] where the regularized solution u is obtained by solving the fourth-order Euler- 
Lagrangian equation. Bredies et al. |15j proposed the total generalized variation regularizer involving a linear 
combination of higher-order derivatives and the TV of u to model the image denoising while Chang et al. 
|25j considered a nonlinear combination of regularizer based on first and second order derivatives. For image 
inpainting, a high order regularization based on Euler’s elastica of u is used in m- Similarly the Euler’s 
elastica energy [56] and mean curvature [MIES] are also proposed to transform the template image T{x + u) 
to map the reference image R{x) in image registration; see also |35||5T]. The above mentioned high order 
regularization methods are effective but due to high nonlinearity efficient numerical solution is a major issue. 

The second type introduces fractional-order derivatives, which are widely studied in other research subjects 
beyond image processing 0 m m m iisi, into regularization of images. For example, Bai and Feng m 
introduced first fractional-order derivative into anisotropic diffusion equations for noise removal 

^ = -D^icilD^^uDD^u) - Z?“*(c(|Z?“n|)i?“w), (4) 

where c(-) denotes the divergence parameter and denotes the adjoint operator of which may be 
viewed as a generalization of the Perona-Malik model. Although the above equation can be related to 
the Euler-Lagrange equations of an energy functional with the fractional derivative of the image intensity, 
generalizing commonly used PDE models, the energy minimization models are not studied as such. The 
discrete Eourier transform is used to implement the numerical algorithm assuming a periodic input image 
at its borders m- See also [45l |44l |4g [66] for more motivations and studies based on the above diffusion 
equation. Chen et al. |28l EH ES] considered the fractional-order TV-L^ image denoising model 

mm {e{u) '■= j^\l + {D^uYdVL + - /II 2 } (5) 

and numerically obtained improved denoising results over the Perona-Malik and ROE models; however no 
analysis was given. There, they converted this primal formulation into a dual problem for the new dual 
variable p = (pi,P 2 ) by m = / — div“p/A and used a dual algorithm using the gradient descent idea similar 
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to the Chambolle method [15] for the ROF. In [5T], the authors proposed a discrete optimization framework 
for image denoising problem where the fractional order derivative is used to model the regularization term, 

N L 

E l(V“<,| + l/2 5]2-2^-*^|[A(/-n),]p, l<a<2,0<s, <l}, (6) 

*j=i i=o 

which is solved by an alternating projection algorithm. See also m- 

These works have reflected good performance of the fractional order derivative in achieving a satisfactory 
compromise such as no stair-casing and in preserving important fine-scale features such as edges and textures. 
These encouraging results motivated us to investigate this new model more closely. 

There have been several other works involving discrete forms of an a-order derivative proposed to tackle 
image registration problem |78[ 154] and image inpainting problem |83| . Comparing with the first type of high 
order models [niisiiin], a fractional order model (type two) is less nonlinear and hence is more amenable to 
developing fast iterative solvers. Clearly there is strong evidence to suggest that fractional order derivatives 
may be effective regularizer for imaging applications. There is an urgent need to establish a rigorous theory 
for the total a-order variation based variational model so that further applications to image inverse problems 
can be considered in a systematic way. 

1.3 Our contributions 

This work is substantially different from previous studies. We mainly focus on the continuous total a 
variation-based model, instead of discrete formulation, and its analysis and associated numerical algorithms. 
Our contributions are four-fold: 

• We analyze properties of the total a-order variation laying foundations for applications to image inverse 
problems as a regulariser; 

• We give a new method for treating non-zero Dirichlet boundary conditions which represents a general¬ 
ization of similar results that existed only in ID to 2D; 

• We establish the convexity, the solvability and a solution theory for the total a-order variation model 
to make it more advantageous to work with than high order and non-convex counterparts (such as a 
mean curvature based model) which are not gradient based and do not have much known theory on 
their solutions; 

• We propose and test four solution algorithms (respectively Split-Bregman based, forward-backward 
algorithm, Nesterov accelerated method and fast iterative shrinkage-thresholding algorithm - FISTA) 
to solve the underlying total a-order variation model. We also compare with related models. 

Our work is hoped to motivate further studies and facilitate future applications of a-order variation based 
regularizer to other imaging problems in the community. 

The rest of the paper is organized as follows. Section 2 reviews the definitions and basic properties of 
the fractional order derivative. Section 3 first defines the total a-order variation and the space of functions 
of fractional-order bounded variations. In this space, it then analyzes the the existence and the uniqueness 
of the solution of the total a-order variation based model for denoising. In Section 4, a boundary con¬ 
dition regularization method for treating nonzero Dirichlet boundary conditions is proposed to effectively 
employ and compute the fractional order derivatives of an image. Section 5 first discusses the discretization 
of the fractional order derivatives by a finite difference method and presents a Split-Bregman scheme for 
effective solution. Section 6 takes the alternative discretise-optimize solution approach and develops three 
optimization-based algorithms (Forward-backward algorithm, Nesterov accelerated method and FISTA) to 
solve the image denoising model. Experimental results are shown in Section 7, and the paper is concluded 
with a summary in Section 8. 


2 Review of fractional-order derivatives 

This section reviews definitions and simple properties of a fractional order derivative which has a long history 
and may be considered as a generalization of the integer order derivatives. Three popular definitions to 
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be reviewed are the Riemann-Liouville (R-L), the Griinwald-Letnikov (G-L) and the Gaputo definitions 

[55115^1^ . 

In this paper, a fraction a € is assumed to lie in between two integers n—l,n i.e. 0 < £ = n— 1 < a < n 
and a fractional a-order differentiation at point a: G K is denoted by the differential operator , where a 
and X are the bounds of the integral over a ID computational domain. Undoubtedly, the gamma function is 
very important for the study of fractional derivative, which is defined by the integral [63] 

poo 

T{z) = / dt. 

Jo 

One of the basic properties is that r(z + 1) = zT^z) and hence r(n) = nl. Before introducing formal 
definitions, we review the following informative but classical example: 

Example 1 The Abel’s integral equation, with, 


1 

fR 



{x 


- t)1-“ 


dr = f{x), 


X > 0 


(7) 


has the solution given by the well-known formula 


tp{x) 


1 d 
r(l — a) dx 



/(t) 



dr, 


cc > 0. 


( 8 ) 


This exarnple helps to understand the formal definitions of fractional derivatives. In fact for 0 < a < 1, 
equation (j^ taking on the form := = f{x) is called the fractional a-order left R-L 

integral of '0(x), and equation (|^ taking on the form ~ is defined as the fractional a-order 

left R-L derivative of f{x). As operators, under suitable conditions [SH], we have x] ~ ^ where I 

denotes the identity operator. 

The first definition of a general order a derivative is the left sided R-L derivative 


D[a,x]fix) 


1 MV 

r(n —a) \dx) (x — ’ 


(9) 


Subsequently the right-sided R-L and the Riesz-R-L (central) fractional derivative are respectively given by 


Dix,b]fix) 


(-1)" MV [' 

r(n —a) \dxJ Jx; (t — x)““"+i 


and 

The second definition is the G-L left-sided derivative denoted by 





fix-jh), 



i(a - 1)... (a - j + 1) 


( 10 ) 


which resembles the definition for an integer order derivative, where [i?] is the integer such that 1 < ["d] < d. 
The third definition is the Gaputo order a derivative defined by 


^'D^a,x]fix) 


1 r 

rin — a)Jx^ (x — t)““”+i 


( 11 ) 


where denotes the n‘^-order derivative of function fix). The right sided derivative and the Riesz-Gaputo 
fractional derivative are similarly defined by 




(-1)" 

r(n — a) (r — x)““"+^ ’ 


^DlaMfi^) = I (^^[a,x]fix) + i-l)^^Dt^,b]fix)) ■ 


When a = n — 1 is an integer, the above left-sided R-L definition reduces to the usual definition for a 
derivative. One notes that when a function is n — 1 times continuously differentiable and its nth derivative is 
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integrable, the fractional derivatives by the above definitions are equivalent subject to homogeneous boundary 
conditions |63| . However we do not require such equivalence for our image function u; refer to Remark]^ 
later. 

Fractional derivatives have many interesting properties — below we review a few that are useful to this 
work. 

Linearity. For a fractional derivative by any of the above three definitions, then one has 

Dia,x]{p f{x) + q g{x)) = p + q 

for any fractional differentiable functions f{x),g{x) and p, g S M. This property will be shortly used to prove 
convexity and to derive the first order optimal conditions. 

Zero fractional derivatives. An integer derivative of an image u at pixels of flat regions may be close 
to zero but the left R-L derivative of a constant intensity function is not zero. One advantage of minimizing 
a R-L derivative instead of the total variation (image gradients) could be a non-constant solution. It would 
be interesting to know the kind of functions that have zero a-order derivatives. 

Lemma 1 (Singularity) Assume that is one of the above three fractional-order derivative operators. 

For any non-integer a > 0 and x > a, there exists a non-constant value function flj) in (a, a;] such that 

Proof. We give explicit constructions. Here we only consider the R-L and Caputo derivatives; for G-L 
derivative, we can derive a similar conclusion through their equivalency. 

1. Assume that 0 < a < 1, for some a; > 0, if /(r) is taken as 

/(r) = {x - 2 t){x - t)°‘ 

for any r G (0, a;] in Abel’s inverse transform then ~ = 0; 

2. Assume that a > 1, if /(t) is taken as 


/(r) = (a; - r)“ ^ 


for any r G (a, x] in a-order R-L derivative, then ~ 

3. Assume that a > 0 in Caputo derivative definition, if /(t) is taken as 


/(r) = (x-r)” 1 


for any t G (a, x] in equation (11), then ^x]f(^) ~ 0- 


Actually for any a > 0, the left R-L x]f(^) = 0 if /(x) = x“ ^ for all fc = 1, 2,..., 1 -|- £ (note £ 
n — 1); refer to [17] . 


[a] = 


Remark 1 For our later applications we take 1 < a < 2. Hence we have the left R-L ^ 

/(x) = x““^ or x““^ i.e. f{x) = x°'® or x~^ '^ when a = 1.6. For the Caputo derivative, ^x]fi^) = 0 */ 
/(x) = 1 or X. 


Boundary conditions. For the left R-L derivative -Dj^^j/(x) of /(x), one assumes that /(a) = 0 or 
f{b) = 0 for the right R-L derivative; otherwise there is a singularity at the end point. So the Riesz R-L 
derivative would require f{a) = f{b) = 0. One solution for nonzero Dirichlet boundary conditions for / 
would be to extract off a linear approximation g(x) (that coincides with f at x = a, b) and to consider 
^[a x]ifi^) ~ 5(^))! however there was no such a method for the 2D case. In Jumarie’s work [50], a simple 
alternative is to modify the R-L derivative to the following 


D“a,x]f(x) 


1 

r(n — a) 



(x - r)“-"+i 


also ensuring that the new fractional derivative of a constant is equal to zero and removing the singularity 
at X = a [8|. In Section]^ we present one method for treating nonzero Dirichlet boundary conditions in 2D. 





Jianping Zhang and Ke Chen 


6 


3 The total o-order variation and its related model 


This section first studies the properties of the total a-order variation, second analyzes a total a-order variation 
based denoising model and hnally presents a numerical algorithm. For the classical total variation based 
model, its solution lies in a suitable space called the function space BV(n) of bounded variation [251[5^fT5] . 
From tests, the total fractional-order variation model can preserve both edges and smoothness of an image; 
we anticipate from the former that its solution should lie in a space similar to the BV space and from the 
latter that the smoothness is due to the non-local nature of the new regulariser. 

It turns out that for total a-order variation using a-order derivatives, a suitable space is the space BV“(n) 
of functions of a-bounded variation on Q which will be defined and studied next. The work of this section is 
motivated by analysis of the total variation (TV) [T1[S1|TH] and of the total generalized variation (TGV) |15) . 

In variational regularization methods, integration by parts involves the space of test functions in addition 
to the main solution space. Before discussing the total a-order variation, we give the following definition: 

Definition 1 (A space of test functions) Let C^(n,K‘^) denote the space of order continuously differ¬ 
entiable functions. Furthermore for any C^(n,M‘^) 9 f i—> if the (i -I- l)th order derivative is 

integrable and ^ [on = 0 for aZH = 0,1,..., £, v is compactly supported continuous-integrable function in 
n. Therefore the l-compactly supported continuous-integrable function space is denoted by 


Definition 2 (Total a-order variation) 


Let K denote the space of special test functions 


K:= |</>G \cf,{x)\ < IforallxGn^ 


where |^| = y ^1- Then the total a-order variation of u is defined by 

TV°‘{u) := sup f ( — u div°‘ 4>]dx, 

^GKJn ^ ' 

where div°‘cf = and denotes a fractional a-order derivative of (pi along Xi direction. 

We note that TV“(m) is the same for any definition of because 4> satishes the equivalence conditions. 
However for our applications in the paper, is generally not the same for different fractional derivatives 
(not even in the distributional sense). 

Based on the a-BV semi-norm, the a-BV norm is defined by 

hllBV° =lklUi+TV“(u), (12) 


and further the space of functions of a-bounded variation on Ll can be defined by 

BV“(H) := {m e L^(H) | TV“(w) <-foo}. (13) 

Lemma 2 (Lower semi-continuity) Let {u^(a;)} be a sequence from BV°‘{Ll) which converge in L^{Ll) to 
a function u{x). Then T\C(u) < liminf T'V°‘{u^). 

k—^oo 

Proof. Since Uk G HV“(H), for any cj){x) G ‘^q(H,IR‘^) such that \(p{x)\ < 1 on H, then div“ cp is bounded, 
hence 


[ ( — u div“ (f>) dx = liminf f 

Jn ^ ' fe-j'+oo 


Uk div“ <f>\dx < liminf TV“(ufe) 

/ fe—>- + oo 


from Uk —> u in L^(H). Taking supj,r^> in the above inequality, we have lower semi-continuity from TV“(rt) < 
liminfTV“(..)(i[Tl[33]forTV?aW 

fe—>-+oo 


Lemma 3 The space 51/“ (H) is a Banach space. 
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Proof. First we can see that BV“(r2) is a normed space following immediately from the definitions of ||M||Li(n) 
and total a-order variation TV“(u), so it only remains to prove completeness. Suppose is a Cauchy 
sequence in BV“(fl); then, by the definition of the norm, it must also be a Cauchy sequence in 
According to the completeness of L^{n), there exists a function u in L^(n) such that —>■ m in L^{n). 

Since {u^} is a Cauchy sequence in BV“(ff), ||u*||bv“ is bounded. Thus TV“(m^) is bounded as fc —)■ oo, 
by the lower semi-continuity of TV“(n) in BV“(f2) space (see Lemma[^, one shows that u € BV“(n). 

We shall show that —>■ u in BV“(n). We know that for any e > 0 there exists a positive integer N 
such that ||m^ — u^\\bv<^(o.) < £ for any j,k > N, hence one has TV“(m^ — u^) < e. Since —?> m in L^(n), 
thus u-l — ^ — u in L^(f2). Hence by Lemma 

TV“(m^ -u)< liminf TV“(m^ -«'=)<£, 

fc—>-oo 

which shows that —>■ it in BV“(f2), therefore BV“(r2) is a Banach space. 

Remark 2 In the literature m, the equivalence of different fractional derivatives requires stringent conti¬ 
nuity conditions e.g. one has b]hi^) = ^6®^ space ‘rfQ~^{[a, 6],IR). However for imaging 

applications (the objective function u), we do not require such equivalence. 


To distinguish the two definitions, we shall continue using the superscript C for C derivatives based 
quantities such as “^div" and ‘^V“ while no superscript means that a quantity is based on the R-L derivative. 

For any positive integer p e N+, let W“(H) = {it e L'Pffl) | ||M||w“(n) < +oo} be a function space 
embedding with the norm 


\u\\w^in) = ( / \u\Pdx -f [ , where V“ii = (^,..., 


9x1 


dxd 


For any ^(x) S Wi{[a,b]) and rj{x) G ^([a,5],I 


/ i{x) ■ ^b\n{x)dx 
J a 

ph 71 i _ j _1 ( ^ 

^(-1)” / ^(x) ■ 

J a 


j=0 


c=b 


(14) 


gives the a-order integration by parts formula (see a)- 

Furthermore, applying (14) twice, we have shown the relationship 


i{x) (^{—4>{x)dx = J 4>{x) ■ \7°‘u{x)da 


(15) 


where it(x) G kF“(H) and (f>{x) G '^g(H,K‘^); clearly the operator (—l)”'^div“ is the adjoint of operator V“. 
Note that, for (f>{x) G ^q(H,IR‘^), we have ^div“0 = div“^ which may not be true if 4>{x) is in a different 
space. 


Proposition 1 Assume that u G 1F“(H), then TV^{u) = |V“it|fix. 


Proof. For any a > 0, using the dual relationship (15), one can obtain that 

f it(x)div“0(x)(ix =(—!)" f (f>{x) ■'V°'u{x)da 

Jq Jn 

and in addition \(p\ < 1 in K implies that 


4>ffx) = 


(-l)”V“w/|V“ii|, 

0 , 


|V“ii(x)|^0; 

otherwise 


can maximize the functional J^<p{x)-V°‘u{x)dx = |V“ii|(ix. By multiplying </)g by a suitable characteristic 
f-compactly supported continuous function in H (e.g., G and then mollifying (see [1] for TV 

and El mi), the new <Ag(x) • V“ii(x)cix with G K is arbitrarily close to |V“M|cix as e —>■ 0 [33], hence 
one shows that TV“(m) = |V“ii|dx by taking sup f^u(x)-div°‘(f^(x)dx = sup (—1)" f>^(x)-V°'u(x)dx. 

4>^^K 4>^^K 
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Remark 3 Since u G W^(i2) leads to TV^{u) = |V“M|fi 2 ;, in fact, it is easy to show that the lower 

semi-eontinuity \W°‘u\dx < liminf \W°‘u^\dx holds in the space similar to the TV ease f,?i9j/ ). 


Lemma 4 The space Wp{fl) is a Banach space. 

Proof. The p = 1 case is clear. Now for p ^ 1, let q satisfy 1/p + 1/q = 1. To obtain the lower semi¬ 
continuity, taking u G Wp{il) and 4’i^) S {0 S ‘^Q(n,R‘^ ) m^)\\L i{n) < 1 for all a; G n}, the following 

inequality 

[ {—l)'^V°‘utljdx= f udw°‘'iljdx = liminf f u^div^^tpdx = liminf [ (— 

Jn an fc-n-oo Jq k-^+oo Jq 

\V‘^u'^\Pdx^ ^ (^J \V°‘u'^\Pdx 


ijp 


holds; further one has |V^uj^da;)< liminf |V“u^|^da;)Then we can deduce the result, follow¬ 
ing the similar lines to proving Lemma 

Lemma 5 The following embedding results hold: lT|'(n) C C BV°‘{C.) C 


Proof. Firstly, from the definitions of BV°‘{n) and Wff{n), we can see that C L^(0) and IFf (fl) C 

L^(n). Secondly, for any / G lT“(n) and cf G K, we have 

f /div“^ dx = (-1)” f cj){x) ■ V°‘f{x)dx < |!V“/||ii(o) < -hoc, 
an an 

i.e., / e BV°‘{n) or ITf(O) C BV°‘{n). Finally W^{VL) C VFf (fl) follows L'^{Vl) C L^{VL). 


Lemma 6 The functional TV^ (u) is convex. 

Proof. The proof follows the linearity of fractional order derivatives, and the positively homogeneous and 
sub-additive properties of TV“ (u). 

Theory for a total a-order variation model. We are now ready to analyze model (HI or the total 
a-order variation model in a more precise form 


mm 

MGBV“(n) 


:=TV“(w) +^F(u)}, F(u) = J \u-z\^dx. 


(16) 


To focus on the total a-order variation model in fl = (0,1) x (0,1) C K^, we assume 1 < a < 2; the following 
theorem establishes convexity of the minimization problem (16). 

Theorem 1 (Convexity) The functional E{u) in BV°‘{n) is convex for A > 0 and strictly convex if X > 0. 

Proof. Since F{u) is a strictly convex functional, the proof follows from Lemma 

If a Banach space X is reflexive (separable), then every bounded sequence in X (in X*) has a weakly 
(weak*) convergent subsequence [501 Prop. 38.2]. Although BV“(fl) is not reflexive, however, it is the dual 
of a separable space. Therefore we can give the following definition: 

Definition 3 (A weak* topology) In BV°‘{Tl), a weak BV^ — w* topology is defined as 


BV^-W 


-G u 


Li(n) 


> u and / 4> ■ V^Uj dx 


(j) ■ V“u dx 


for all (p in (fl, M^^). 

From the above definition]^ we may derive the weak compactness of BV“(n) on the weak* topology. 
This, combined with the weak lower semi-continuity of E{u) and boundedness of Banach space BV“(11) (i.e, 
u is bounded in Banach space BV“(n)), yields the following result: 

Theorem 2 (Existence) The functional E{u) : BV'^{Tl) —>■ K has a minimum. 

Proof. Follow the similar lines of [501 Prop. 38.12(d)]). 

Theorem 3 (Uniqueness) The functional E{u) has a unique minimizer in BV*{Vl) when A > 0. 

Proof. The convexity result of Theorem leads to uniqueness of solutions. Refer to (501 Theorem 47C]. 

We remark that similar existence and uniqueness theories of the total variation problem can be found in 

[Ulllllg]. 
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4 Nonzero Dirichlet boundary conditions and regularization 

The standard definitions for fractional derivatives require a function to have zero Dirichlet boundary condi¬ 
tions due to end singularity, but for imaging applications such conditions are unrealistic and too restrictive. 
To obtain the system for finding the unknown intensities u at inner nodes of a discretization grids, we have 
to use boundary conditions, but the difficulties caused by them in fractional derivative computations would 
be hard to overemphasize; inaccurate boundary conditions can easily lead to the oscillations near boundaries, 
so proper treatment of the boundary conditions for problems involving fractional derivatives is crucial. 

In this section, we shall reduce nonzero Dirichlet boundary conditions to zero ones so that standard 
definitions and our introduced algorithms become applicable. The basic idea of boundary regularization is 
to introduce an auxiliary unknown function which satisfies the zero boundary conditions. In this way, the 
non-zero boundary conditions move to the right-hand side of the equation as a new known quantity. 

We recall that, in the ID case, if the boundary conditions are nonzero 

m( 0) = a, m( 1) = b, 

we can reduce them to zero boundary conditions by introducing an auxiliary function e{x) = a(l — cc) -I- bx. 
Precisely taking u{x) = u{x) — e{x) |Mj, then 

u(0) = 0, u(l) = 0; u'iO) = u'{l) = 0 

and a Neumann boundary condition is imposed by artificially extending the boundary values i.e. e'(0) = 
e'(l) = 0 on dfl. 

Below we generalize the above ID idea to the 2D case, assuming that the four corners of the solution are 
given or accurately estimated: 


u(0, 0) = a, u{0, 1) = 6 , m( 1, 0) = c, u(l, 1) = d. 

With a,b,c,d known, at any image point {x,y) € D, a bilinear auxiliary function satisfying the above 4 
conditions ei{x, y) = a + {c— a)x -I- (6 — a)y + {d + a — c — b)xy can be constructed to lead to 

u{x,y)=u{x,y)-ei{x,y) (17) 


which takes zero values at all 4 corners. 

If boundary conditions u(0,y) = ai{y), u{l,y) = 02(2/), w(a:, 0) = bi(x), u{x, 1) = 62(2:) at dC. are known 
a priori, then we can easily verify that 

di{y) : = u{Q,y) = ai{y) - ei(0,2/), 02(2/) := u{l,y) = 02(2/) - ei(l,2/); 

bi{x) : = u{x, 0) = 6i(a;) — ei(a;, 0), b 2 {x) := u{x, 1) = 62(2;) — ei{x, 1) 

define the new Dirichlet conditions for u{x,y). 

We can achieve zero conditions at the edges using the auxiliary function 62(2;,2/) = ^(1 — x)di{y) + 

2^02(2/)^ + ^(1 — y)bi{x) -I- 2/62(2:)^. It is clear to see that the new image u{x, y) = u{x, y) — ei(x, y) — e 2 {x, y) 
satisfies 

u{x,y)\dn = 0- ( 18 ) 

Remark 4 It remains to address the question of how to obtain estimates ofu{x,y) at corners and edges: 

1. The true intensities a := u(0,0), b := m(1,0); c := tt(l,0), d := u(l,l) in four corner points are 

unknown a priori, to build the auxiliary function ei{x,y), the solutions approximating to them should 
be solved from the observed image z{x,y) by the local smoothing or other simple techniques. 

2. Similarly, the true edge intensities ai{y), 02(2/), ^1(2:) and 62(2:) are also not given, a reconstruction 
step on boundary dfl must be proceeded in order to capture a robust solution. To do this, we can apply 
a ID model. 

According to Remark we can propose a complete procedure for regularizing boundary conditions for 
2D variational image inverse problems in edges and corners. 
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• Firstly, we restore iniage intensities in 4 corner points from an observed image z. A natural technique 
would be local smoothing operator for the region of corner points, the oscillations could also be reduced 
by many variational methods to local regions. 


Secondly, in order to reconstructed 4 edges from the restored intensities z{0,y), z{l,y), z{x,0) and 
z(x, 1), the total a-order variation regularization would be used to solve four ID inverse problems i.e. 
solve an equation like (16): 



(19) 


Thus through ei{x,y),e 2 {x,y), we see that equation 0 reduces to finding the new image u{x, y) with zero 
Dirichlet conditions and hence the standard definitions of fractional derivatives for u(a;, y) apply. 


5 Discretization and Split-Bregman algorithm 

Since solution uniqueness of our variational model 0 is resolved, we now consider how to seek a numerical 
solution of the total a-order variation model. We first reformulate it in preparation for employment of 
an efficient solver and then discuss some discretization details (by finite-differences) before presenting our 
Algorithm 1. 


5.1 A Split-Bregman formulation 

Inspired by Goldstein and Osher’s Split-Bregman work I43|. we introduce a special and new variable d{x) = 
{di{x), d 2 {x))’^ to the total a-order variation based model (161 to derive the following constrained optimization 
problem: 


min f \d\dx+^F{u), s.t. d = V“u. 

u,d Jq 2 

To enforce the constraint condition, we transfer it into the Bregman formulation 

, d'^'^^) = min f \d\dx+^F{u)— f <p^,d—d^>dx 
“.d Jo. 2 Jq 

dx + ^ f |d — 

2 Jn 


( 20 ) 


- <Pu-,U-U^ > 




The above iterative scheme can be simplified to the two-step algorithm [la [701 In]: 


min [ \d\dx + ^ [ \d — V°‘u + —\'^dx + [ \p\'^dx + ^ F(u) 

“A Jn 2 Jn h Jn 2 


( 21 ) 


with the multiplier updated by iteration = p^ — y(d — V“u), where p{x) = {pi[x),p 2 {x))'^■ Here the 
two main subproblems of (21) are 


|d — -I- —pdx; 


Subproblem d : min / \d\dx + ^ 

d- Jn 2 Jn 

Subproblem u : min J(u) := ^ / \d — V°‘u + —\'^dx + ^F{u). 

u 2 Jn p 2 


( 22 ) 


Further note that the subproblem d has a closed-form solution [48] . while the subproblem u is determined 
by the associated Euler-Lagrange equation as shown below. 


Theorem 4 Let u{x) he a minimizer of functional J{u) from (22). Then u{x) satisfies the following first 
order optimal condition 

\ ^ 


+ A(w — z) = 0 


(23) 
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with one of these sets of boundary conditions 

du{x) 


i) fixed : = &i(x), and 


dn 


an 


= hix); 


ii) homogeneous : D°‘ ^ [ V“u — d — — ) • n 

h 


P 


an 


= 0, 1 -n 


= 0 


an 


where n denotes the unit outward normal and div^ denotes the divergence operator based on the C derivative. 
Proof. Refer to Appendix. 

5.2 Discretization of the fractional derivative 

Before introducing the finite difference discretization of the fractional derivative, we define a spatial partition 
{xk: ( for all fc = 0,1,..., + 1; Z = 0,1,..., M + 1) of image domain Assume u has a zero Dirichlet 

boundary condition (practically we apply the regularization method ^ first before discretization). Here 
we mainly consider the discretization of the a-order fractional derivative at the inner point {xk,yi) (for all 
fc =, 1,..., A^; Z = 0,1,..., M) on ft along x-direction by using the approach 


D[a,b\fixk,yi) = 


^Uixk,yi) 


h° 


+ Oik) =-( 


1 /'S‘ff{xk,yi) , S‘ff{xk,yi) 


k+i .. . 

3=0 ■ “ 


h° 

N-k+2 


h° 


Oih) 


(24) 


3=0 


which is applicable to both the R-L and C derivatives 
j = 0,1,..., A^ + 1 and 


Eg, where fi = fs,u = (-1) 


w, 


(“) 


= iCd =(1 - 


1 + a, 


)w|“\, for j > 0. 


Alterative discretization for fractional derivatives in the Fourier space can be found in [13 m]. 

Observe from (24) that the first order estimate of the a-order fractional i,]fi^k,yi) along x-direction 


at the point {xk, yi) with a fixed yi is a linear combination of A^ -I- 2 values {/g j /li • ■ • I In^ fN+i}- 

After incorporating zero boundary condition in the matrix approximation of fractional derivative, all N 


equations of fractional derivatives along x direction in (24) can be written simultaneously in the matrix form 
(denote w = ujq + uj^)- 

( fi \ 

n 


/ Sgf{xi,yi) \ 


/ 2a;? 

W ujf ■ ■ ■ U3fq 

\ 

Sof{x2,yi) 


VJ 




1 





“ 2Zi“ 


. .a 
Wg 


V SUixN,yi) ) 



2a;? w 



V 

uj^ w 2a;? 

/ 


(25) 


V /v / 


From the definition of fractional order derivative (24), for any 1 < a < 2, the coefficients u!^'^ has the 
following properties [63l ES] : 


1). 

OJq 

= 1, a;g“^ = —a < 0, 

2). 

1 > 

3). 


o' 

II 

3 

o 

4). 

Ek=o 


> > 
(“) ^ n 


> 0 , 


Hence by the Gerschgorin circle theorem, one can derive that matrix in (25) is a symmetric and negative 
definite Toeplitz matrix (i.e. is a positive definite Toeplitz matrix). 

We recall that the Kronecker product AC) B oi the px q matrix A = [a^] and the nx m matrix B = [brt] 
is the np x mq matrix having the block structure A® B ■.= [aijB]. Further vector (A® B)x can be computed 
by matrix scheme BXA^ (i.e., [(A ® B)x]s = [BXA^]j^i with s = {i — l)m + j), where the m x q matrix X 
is the reshape of the vector x along its column. 
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Let U G denote the solution matrix at all nodes {khx'Jhy), k = 1,N; I = cor¬ 
responding to x-direction and y-direction spatial discretization nodes. Denote hy u G the ordered 

solution vector of U. The direct and discrete analogue of differentiation of arbitrary a order derivative is 

where ■ ■ ■, ..., , u = {un,..., mati, wi2, • ■ ■, unm)^ ■ Similarly, the a-th order 

y-direction derivative of u{x]y) is approximated by: 


4“^ = 0 Jjv)u, where 


(u 


(a) 

11 I•• ■ I 


^IM’ 


(a) (a) \ 

M21 , ■ ■ ■ , J 


T 


5.3 The Split-Bregman algorithm 

In discrete form, we are ready to state the discretized equations in structured matrix form. The discrete 
scheme of (231 is given by 

■ =0 


1 

h 


with discretizations di = (d^, ■ ■ •, ^12: ■ • ■: and p, = {p\^,... ,pVnPi2. ■ • • ^PnmV vectors d 

and p (i = 1,2). A matrix approximation equation is given as 


iB%U) + U{B^)^ B‘^]+XU = XZ+[ {B%)^ D, + ^2^^ + - (B^)^ Bi + B2B 




-“M ) I 


(26) 


wu 


where Di and Pi are N x M-size reshape matrices of vectors di and pi for i = 
following justifies the use of a conjugate gradient method for WU = F. 

Theorem 5 The weighted matrices inner product {WU,U) = J2iJ2WikUkj)Uij 

ij k 

U ^ 0, where W is a known positive definite operator. 

Proof. For any matrix B 7^ 0, it is easy to show that 


1,2, A = (-l)"A//r. The 
is positive for any matrix 


{WU,U) = (((B^)^(B^t/) + C/(B^)^B^) +AC/, u) 

= {B%U, B%U) + (C/(B^)^, t/(B^)^) + X{U, U) 
= \\B%U\\l + \\U{Blj)'^\\l + X\\U\\l>D, 


which completes the proof. 

An implementation of this method may be summarized below: 
Algorithm 1 (Split-Bregman iterations (PDE-SB)) 

step 1. Boundary regularization for an observed image z; 
step 2. Given initial matrices Pi^^, B ^and 


step 3. Solve subproblem d.' Compute the auxiliary matrix 


Di 

D 2 


from the closed form solution 



k+l 


= shrink 


B“C/'=+i 



by solving the Moreau-Yosida problem with the f regularization; 
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step 4- 
step 5. 


Solve subproblem u: Find the solution of (26) with an effective parameter \ p. by CG method; 

fc+1 / \ k // \/\ ^+1^ 


Update 


Pi 

P2 


Pi 

P2 


■7 




Di 

D2 


with 7 G (0,1]; 


step 6. Check the stopping condition; 

• If\U^ - U^+^\ < e, 

stop and return U* := C/^+^; 

• else 

k := k + 1, go back to 

• end 


step 1. Accept the correct solution U from boundary regularization. 


6 Optimization based numerical methods 

As many variational models are increasingly solved by the discretise-optimise approach, we now present three 
related algorithms for model Q after applying a finite difference discretization. In this section, we assume 
that we have the zero Dirichlet boundary conditions for u mainly to simplify the notation. 

As in ^ the a-th order derivative of u{x; y) along all x-direction nodes in kl can be given by matrix 
BfjU, and similarly U{B’fj)'^ for j/-direction (as U is the solution matrix). 

Define {U,V) = YlUijVij and let lA = {p | 0 < p < 1}, lA = {p | |p| < 1}. Then using the discrete 
ij 

setting introduced above, the discretised problem of model ( [l^ is 

uiinmaxG{U,D*<h) + ^H{U) (27) 

ueVi $gV2 2 

where H{U) = and G(t7,D*$) = ([/,D*$) = {Pn^i P due to = 

B%^i -h ^ 2 {B‘fj)'^. We also have the adjoint relationship {U,D*^) = {DU,^) with DU = {BffU,U{Bfj)'^) 
and $ = (<i>i, <i>2). In line with the literature, this model can be denoted by the convex optimization problem 
in a generic notation by 


min{/i(x) +/2(x)} i.e. min{/i(x) +/2(y)} s.t. x = p 

X x,y 


(28) 


where one views x = U, /i(x) = max<jgV2 G{U, 77*$), /2(x) = H{U). We also need the notation 


prox)^(x'=) := argmin|/i(x) + | 


where fi can_te any other convex function and A > 0. 

To solve (28) by the methods to be presented, computation of the proximal point proxj^(x) is a major 
and nontrivial step. We consider how to compute it when D = V“, borrowing ideas from solving a similar 
problem of TV regularization. In a dual setting, Chambolle [MIID] firstly proposed a discrete dual method 
by optimizing a cost function consisting of two variants [niiiH]. Recently, one variant of this scheme is 
employed in |28j to effectively solve a fractional image model by a dual transform. The other variant is used 
in [26]. 

Define two projections as 


Projy Ip) = 


0 

P 

1 


p < 0 
0 < p < 1 

1 < P, 


Projy^b) = 


max(l,||p||)' 


Noting |-]^g optimal solution is 


X = prox]^(x'') = Projy^(x) 


(29) 
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where x = and d* is unknown. Based on methods of Chambolle m and Beck-Teboulle m , we 

see that (29) can be used to reduce the min-max problem 


min max(x, llx — x 


k\\2 


to the dual problem max^^y^(x), + ^||Projy^(x) — x^|p and further to 


(Projy^(a;),D*$) + — ||Projy^ (x) - a;'^f = (Projy^ (x), £>*$) 

+ ^||Proj,.^(i) -x^ + - ^\hD*n^ - ^2(Proj^^(x) - x^-fD*^) 

= ^l!Projy,(x) - {x’^ - 7D*^)r - ^\\lD*n^ + ^2{x\^D*<^) 

= ^|!Proj^^(x) - (0.'= - ^ {\\lD*n^ - 2{x^^D*^) + \\x^r) + 

= ^\\PT0iy^{x) - {x^ - - ^\\x^ - 7i^*$f + ^Wx'^r 


1 


|x-Projy^(x)f - ||xf+ ||x 


k\\2\ 


(30) 


i.e. max$gV2(Projy^(a;),£>*$) + ^||Projy^(a;)-a;'=|p = min /i(d>) where h($) = \\x’^ - 'yD*^\\^ - \\{x'^ - 


4>eV2 


711*$) - Projy^(a:'= - -/D*<h)\\‘^ - \\x'^p = ||a;|p - ||a: - Projy^(a;)|p - ||a:'=p. 

Below we consider the operator S{x) = ||x — Projy^ (^)ll^ = inf |<5 v^l( y) + ^lly ~ Since its gradient 

is \7xS{x) = 2(x — Projvi{x)), we get 


V$h($) = -27i?(Proj^^(x'= - 721*$)). 


The minimization problem min /i($) can be solved to obtain the $-update as follows 

•s>eV2 


1) $ = - L(/i)V$/i($”); 

2) $"+i = Proj^^($) = ProJv, ($” + 2L(h)7Zl(Proj^^ - 7D*$"))) 


using the gradient projection scheme of /i($) [13]. Here L(h) < 167^ is the Lipschitz constant. Finally the 
proximal point prox^^(a;^) is given by (29) once $ is obtained; see also [T3] . 


6.1 Forward-backward algorithm 

Various applications in sparse optimizations stimulated the search for simple and efficient first-order methods. 


The forward backward scheme for (28) is based (as the name suggests) on recursive application of an explicit 
forward step with respect to /2, i.e, 


f2{x’") + {Vf2{x^),x- : 

-V- 

l(x) 




27 


and followed by an implicit backward step with respect to /i, i.e., 

mm{/i(x)-b }. (31) 

The scheme decouples the contributions of the functions fi and /2 in a gradient descent step [12] . The scheme 
is also known under the name of proximal gradient methods [TusnuToidi], since the implicit step relies on 
the computation of the so-called proximity operator. 

The forward backward algorithm is summarised as follows. 

Algorithm 2 (Forward-backward algorithm IFBI [12] 1 
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• Fix initial Xq, set e G [0,min{l, 1//?}], /3 fa Lipschitz parameter); 

• For k > 0 

Step 1. jkG [e, 2/^ - e], \k G [e, 1]; 

Step 2. Uk = prox'l’°{xk) 

Step 3. Xk+i = prox'jfiyk); 

Step Xk-^-l Xk “1“ ^k(^Xk-i^l 3^/c); 

Step 5. Stop when ||xfc+i — Xk\\ is small enough otherwise continue. 


6.2 Nesterov’s method 

As a gradient based method, though simple, the above method can exhibit a slow speed of convergence. 
For this reason, Nesterov m proposed an improved gradient method aiming to accelerate and modify the 
classical forward-backward splitting algorithm, while achieving an almost optimal convergence rate. As a 
consequence of this breakthrough, a few recent works have followed up the idea and improved techniques for 
some specific problems in signal or image processing [niin]. 

Recently Nesterov [5S] presented an accelerated multistep version, which converges as 0(i) (r is the 


iteration number). For a problem of type (28), this new method introduced a composite gradient mapping. 
We now show the algorithm as follows. 

Algorithm 3 (Nesterov accelerated method (Nesterov |58| 11 

• Fix initial Xq, bo, set j/g = xq and P (a Lipschitz parameter); 

• For k > 0 

Step 1. Find a = Ok from the quadratic equation 2 (i^,,+a) ~ ’ 

Step 2. v=prox^){xk - Pk); 

Step 3. Zk+i = '’"bt+af" 

Step 4 . Xk+i = proiPf^ {zk+i - f 2 [zk+i))i 

Step 5. pk+i =yk + afeV/2(a:fe+i); 

Step 6. hk+i =bk + ttfe/ 

Step 7. Stop when ||a;fc+i — Xk\\ is small enough otherwise continue. 


6.3 FISTA method 


Beck and Teboulle mill] proposed a fast iterative shrinkage thresholding algorithm (FISTA) to solve the 
image denoising and deblurring model, The method applies the idea of Nesterov to the forward-backward 
splitting framework, resulting in the same optimal convergence rate as Nesterovs method but wider applica¬ 
bility. It can be applied to a variety of practical problems arising from sparse signal recovery, image processing 
and machine learning and hence has become a standard algorithm. 

Applying it to ( [^ , we obtain Algorithm]^ below. 

Algorithm 4 (FISTA (Beck-Teboulle |12L I13L I14| II 

• Fix initial xo, set zq = xq and to = 1, j3 (a Lipschitz parameter); 

• For k >0 

Step 1. yk= Zk- /3“^V/2(z/c) 

Step 2. Xk+i = proxPj^ (yk) 

oj. o + l + v^ 

Step 3. tk+i = —— 

Step 4 . Zk+i =Xk + {l + ^^){xk+i - Xk); 

Step 5. Stop when ||xfc+i — Xk\\ is small enough otherwise continue. 
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7 Numerical results 


Finally, we present some numerical results from using the four presented algorithms denoted by 


PDE-SB: 

Opti-FB: 

Opti-Nesterov: 

Opti-FISTA: 


PDF-based Split-Bregman (Algorithm [^ ; 

Optimization based Forward-backward (Algorithm]^; 
Optimization based Nesterov Accelerated method (Algorithm]^; 
Optimization based FISTA (Algorithm]^, 


and their comparisons with related methods. In all tests, an initial solution is the noisy image z(x,y), 
the algorithms solving the diffusion equation or optimization problem are stopped after achieving a relative 
residual of 10“^ or a relative error of 10“® within 1000 outer and 15 inner iterations. Here we mainly compare 
the solution’s visual quality, the snr (the signal-to-noise ratio) and psnr (the peak signal-to-noise ratio) values 
which are given 


snr{u,u*) = lOlogj^o 


\u — meamu 


UxTl. 


\U — U’ 


\l 


; psnr{u,u*) = lOlogj^Q 


,(inaxu*^^ 


\\u — W 




where mean{u*) is an average value of the true image u*, Ux and Uy denote the size of the test image z. It 
should be noted however that these valuations do not always correlate with human perception. In real life 
situations, the two measures are also not possible because the true image is not known. 

In general, an optimization problem may be solved many times to select a suitable regularization parameter 
A or to optimize the solution for the underlying inverse problem; a solution is accepted when some stopping 
criterion is satisfied. It remains to carry out a systematic study on our new model as in [82] for the TV 
model. However we shall use the best (numerical) A for all models in the following tests. 

For denoising, F{u) = (u — z)^ is the measure between the solution u and the observed image z. To 
intuitively describe the denoising ability, four sets of data will be used in this part {also see Fig§: 

PI: Problem 1 - Parabolic surfaces; P2 : Problem 2 - Saddle surface; 


P3: Problem 3 - Pepper; 


P4 : Problem 4 - Penguin. 



Figure I: Test datasets. 

Though our framework is readily applicable to image deblurring and image registration, here we only present 
denoising results. 

7.1 Performance comparisons of boundary regularization 

We first test the idea from One the hand, the variational framework seeks the boundary conditions of a 
nonzero Dirichlet or a Neumann type on dfl and also real images do have nonzero boundary conditions. On 
the other hand, fractional order derivatives require homogeneous boundary conditions (as used in works of 
many authors) due to end singularity. In order to aid accurate computation of the discretized fractional order 
derivative, in our work, a boundary processing technique Qhas been proposed to transform nonzero boundary 
conditions of observed data z into zero boundary conditions; hence a consequent matrix approximation to 
the fractional derivative operator can use a zero Dirichlet boundary condition. 

Here we test the performance and effectiveness of our boundary regularization against no regularization. 
The experiment is carried out on PI - Parabolic surfaces as shown in Fig. i.e., a synthetic image of size 
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256 X 256 and range [0, 1], in Fig. [73a) , which is added zero mean value Gaussian random noise with a 
mean variance (5 = ^ to get the noisy image displayed in Fig. |7.lfd). For the boundary regularization 
case, the approximation u|ao from the observed data z\aQ, is from applying one dimensional fractional order 
variation model as described in 0 The treated case is named as ‘Treated’ whose results are depicted on 
Fig. |7.1[ b) and Fig. |2(b)[ ), where psnr= 47.53 and snr=35.43. The solution obtained from assuming zero 
boundary conditions for u is named as ‘Non-treated’ with its results depicted in Fig. (He) and Fig. |2(a)[ 
where psnr= 23.69 and snr= 10.38. Clearly our boundary regularization treatment is effective. 





(a) Slice for non-treated - Bad. 


(b) Slice for treated - Good. 


Figure 2: Test for PI —Comparisons between treated and Non-treated cases for non-zero boundary conditions 
{S = using PDE-SB. The treated case has psnr=47.53 and snr=35.43, while the non-treated case has 
psnr=23.69 and snr=10.38. Clearly our boundary regularization Qis effective while direct application of a 
fractional model leads to incorrect boundary restoration. Here the error r = ||m — u*||j’/||u*||i^’. 


7.2 Comparisons of Algorithms 1—4 

In Table we compare the restoration quality (via psnr and snr) of 4 Algorithms. There, all four test 
datasets are used. In the cases of synthetic images PI and P2 with noise variation <5 = A is taken as 
12000 and 3800 respectively and a = 1.6. In the cases of natural images P3 and P4 with noise variation 
S = A is taken as 18000 and 20000 respectively and a = 1.4. One can see that, from Table Opti- 
Nesterov and PDE-SB perform similarly in terms of the best restoration quality (via psnr and snr). However 
in efficiency (computation times cpu(s)), Opti-FISTA and PDE-SB are the best while Opti-Nesterov takes 
more computational times than other three algorithms. Evidently, overall, PDE-SB (Algorithm 1) shows the 
most consistence in good performance in tested cases considered. 
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Table 1: Comparisons of optimizing Algorithms, where ^ ^ for saddle and parabolic surfaces, <5 = ^ for 

pepper and penguin images. 




Opti-FB 


Opti-Nesterov 

Opti-FISTA 


PDE-SB 



snr 

psnr 

cpu(s) 

snr 

psnr 

cpu(s) 

snr 

psnr 

cpu(s) 

snr 

psnr 

cpu(s) 

PI 

36.78 

50.09 

16.83 

36.91 

50.22 

27.23 

36.94 

50.24 

16.71 

36.96 

50.27 

14.53 

P2 

31.04 

53.08 

17.28 

31.61 

53.67 

28.43 

31.46 

53.50 

18.14 

31.63 

53.69 

15.09 

P3 

29.21 

43.29 

15.96 

29.40 

43.49 

16.09 

29.40 

43.49 

9.75 

29.48 

43.56 

8.16 

P4 

25.19 

38.01 

14.68 

25.35 

38.15 

16.27 

25.34 

38.15 

8.62 

25.34 

38.14 

8.45 


7.3 Sensitivity tests for a and A 


Since our model ( |16[ ) contains two main parameters: a for the order of differentiation and A as the coupling 
parameter for a regularized inverse problem, it is of interest to test their sensitivity on the restoration quality. 
Here we shall test all Algorithms’s sensitivity using the image P2 - Saddle surface of size 256 x 256, after 
adding zero mean value Gaussian random noise image of range [0, 1] and 6 = 

Varying A in a large range from 400 to 60000, all four algorithms are tested on this synthetic image 
with the results shown in Figs. 3(a) and 3(c) for different stopping criterions (GSC: the general stopping 
criterions with the relative residual 10“^, relative error 10“®, inner iterations 10,SSC: the strong stopping 
criterions with the relative residual 10“^, relative error 10“^°, inner iterations 25 ). Different from the TV 
denoising case where the regularization parameter A is crucial for restoration quality [82) . however. Figs. 3(a) 
and |3(c) show that our total a-order variation regularization model still obtains a satisfactory solution for a 
large range of A; this is a pleasing observation. Of course there exists an issue of an optimal choice. 

Next varying a S (1,2) from 1.1 to 1.9, Figs. |3(b)| and |3(d)| show four algorithms’s restored results 
responding to two stopping conditions GSC and SSC. As represented, the smaller a leads to the blocky 
(staircase) effects in u and the larger a will make solution u too smooth along xi- and a; 2 -directions respec¬ 
tively. For denoising, our test suggests that a = 1.6 is suitable for smooth problems because the diffusion 
coefficients are almost isotropic in all regions, leading to smooth deformation fields, and a = 1.4 is appropri¬ 
ate for nonsmooth problems because the diffusion coefficients are close to zero in regions representing large 
gradients of the fields, allowing discontinuities at those regions. 

We should emphasize that the stopping criterions have impacted on the actual numerical implementation. 
In other words, if we drop the limit on the maximal number of inner iterations and relative residuals (and 
relative errors), some methods take too long but obtain the more satisfactory results. 



(a) GSC'.psnr vs. A 


(b) GSC:psnr vs. a 


(c) SSC:psnr vs. A 



(d) SSCipsnr vs. a 


Figure 3: Sensitivity test of Algorithms 1-4 to parameters A (with the fixed a = 1.6) and a (with the fixed 
A = 3800) in the cases of the GSC and SSC stopping conditions. 


7.4 Comparisons with other non-fractional variational models 

In this test, we compare our total a-variation model(PDE-SB) with three popular methods for variational 
image denoising. The first compared approach is naturally the TV model proposed by Rudin et al. m 
because the total a-order variation model in this work is inspired by it. The second compared work is the 
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mean curvature model m which also addresses the problem of restoring a good result for a smooth image; 
their approach is different from ours since it is focused on higher order regularization and a multigrid method. 
See also HSlIITlIll]. The third compared approach is the TGV model m involving a combination of first 
order and higher-order derivatives to reduce the staircasing effect of the bounded variation functional. 

In Table we first compare the restoration quality (via psnr, snr) and efficiency (computation times 
cpu(s)) of four approaches by testing the artificial images (PI - Parabolic surface, P2 - saddle surface) and 
the natural images (P3 - Pepper, P4 - Penguin); in each approach relevant parameters are shown in Table 
We see that, with the emperically optimal parameters A*, the differences of four models are very small, 
though our new and convex model is slightly better. In other tests where such optimal parameters are not 
used, our new model performs more robustly and better. 

In order to present more visual differences, some stronger regularization parameters (A*/2) and higher 
noise variations (with the noise level 5 = ^) are tested, the solution’s visual representations restoring the 
natural image P3 - Pepper in Fig. |4(b)| are shown in Fig. While ROF denoising leads to blocky results, the 
mean curvature model performs better in the smooth regions but exhibits more smooth near discontinuities, 
the total generalized variation model leads to further improvements over the aforementioned models. The 
total fractional-order variation model leads to significantly better results. The reason is that the new model 
tries to approximate the image based on affine functions or non-local high order smooth functions, which is 
clearly better in this case, in other words, our approach is more effective in eliminating the noise for smooth 
images and is competitive to high order methods; in efficiency the new approach (PDE-SB) is much faster 
than the TGV and the mean curvature. We also plot four error results between the restored and true images 
along a diagonal (magenta) line in Fig. |4(a)| for comparison in Fig. we see that PDE-SB produces the best 
restored surface, which show a major advantage (or better performance) of using our total a-order variation 
model (161 when the test image is smooth, and even when the contrast between meaningful objects and the 
background is low. 


Table 2: Gomparisons of four models in restoration quality: the total a variation model (16), Mean Gurvature 
[S31 [T71 [75] , TV [5711311 HSj and TGV [T5| models for synthetic images (PI and P2 in Fig. 0 and natural 

images (P3 in Fig. 
p = 1.1, 7 = 1, a 




with different noise variances 5j = ^ (correspondingly we use \j). We first fix 
^T 6 for P1-P2, a = 1.15 and a = 1.1 for P3-P4 respectively in the total a-order 
variation model, and 7 = 19, /? = 10“® in the mean curvature model, and two weight parameters of the 


first and second order term in TGV model (i^o = 1, ni = 2 for P1-P2, i^q = 1, = 0.5 for P3-P4), other 

parameters are as shown on the “para” rows. A^^ from (19) is required by the new model only. 




Mean curvature 1 75 II 

TV Eg 

TGV [Tg 

Total Q-order model |16)| 


5 

snr psnr 

snr 

psnr 

snr 

psnr 

snr psnr 


10 

33.44 46.74 

32.17 

45.52 

36.41 

49.72 

37.55 50.86 

PI 

20 

30.19 43.50 

29.55 

42.83 

33.03 

46.33 

33.52 46.83 



Ai = 1/0.4 X 256^ 

Ai = 

1026 

Ai = 

1/1.2 X 256^ 

aJ° = 0.1, Ai = 21900 


Para 

A2 = 1/0.03 X 256“’ 

A2 — 

535 

A2 — 

1/0.6 X 256““ 

A“° = 0.1, A2 = 14400 


10 

27.27 49.31 

23.09 

45.13 

30.75 

52.68 

32.02 54.18 

P2 

20 

22.88 44.92 

19.45 

41.49 

25.62 

47.51 

26.48 48.54 



Ai = 1/0.9 X 256““ 

Ai = 

883 

Ai = 

1/0.9 X 256““ 

a;““ = 1 , Ai = 1800 


Para 

A2 = 1/0.01 X 256““ 

A2 — 

488 

A2 — 

1/0.5 X 256^ 

XY" = 0.2, A2 = 1800 


10 

20.43 38.80 

20.08 

38.35 

20.40 

38.78 

20.48 38.86 


15 

18.76 37.11 

18.01 

36.69 

18.68 

37.12 

18.84 37.20 

P3 

20 

17.48 35.82 

17.17 

35.33 

17.55 

35.87 

17.57 35.90 



Ai = 1/16 X 256““ 

Ai = 

2216 

Ai = 

1/55 X 256“ 

a;““ = 1 , Ai = 16500 


Para 

A2 = 1/14 X 256^ 

A2 = 

1373 

A2 = 

1/26 X 256“ 

XY = 0.1. A2 = 9300 



A3 = 1/6 X 256^ 

A3 — 

893 

A3 — 

1/12 X 256“ 

XY = 0.01. As = 6200 


5 

25.16 37.95 

24.85 

37.58 

25.39 

38.20 

25.34 38.14 


10 

21.72 34.60 

21.33 

34.07 

21.82 

34.71 

21.75 34.62 

P4 

15 

19.26 32.05 

18.66 

31.29 

19.44 

32.21 

19.42 32.20 



Ai = 1/9 X 256““ 

Ai = 

3341 

Ai — 

1/49 X 256“ 

aJ° = 0.1, Ai = 24000 


Para 

A2 = 1/5 X 256““ 

A2 = 

1856 

A2 = 

1/20 X 256“ 

xY = 0.1, A2 = 8000 



A3 = 1/6 X 256““ 

A3 = 

1095 

A3 = 

1/11 X 256“ 

xY = 0.1, As = 18500 


8 Conclusions 


The total a-order variation regularization with fractional order derivative is potentially useful in modeling all 
imaging problems. In this paper we analyzed rigorously a simple variational model using total a-order varia- 
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(a) True (b) Noise /. (c) TV. 






(d) Mean Curvature. 


(e) TGV. 


(f) Our Approach. 


Figure 4: Comparison I —Comparisons of our PDE-SB with TV, mean curvature and TGV models. 


tion for image denoising. One Split-Bregman based algorithm and three optimization-based algorithms were 
developed to solve the resulting image inverse problem. Instead of using the usual fixed and zero boundary 
conditions, we proposed a boundary regularization method to treat the fractional order derivatives. Numeri¬ 
cal results show that the PDE-based Split-Bregman algorithm (PDE-SB) performs similarly to (though more 
stably than) optimization-based approaches while our boundary regularization method is essential for getting 
good results for imaging denoising. Moreover, PDE-SB outperforms currently competitive variational models 
in terms of restoration quality. There are still outstanding issues with our proposed model and algorithms; 
among others optimal selection of A is to be addressed. Future work will also consider generalization of this 
work to other image inverse problems. 


Appendix: Proof of Theorem 

To shorten the proof, let w be a function in IV“(D) to be specified shortly. For u G IV“(D) C BV“(D), we 
compute the first-order G-derivative (Gateaux) of the functional J(u) in the direction w by 


,, . J(u-\-tw) — J(u) Qiu-\-tw) — Q{u) \ F(u-\-tw) — Flu) 

J [ujuj = lim -- = hm -h — -—■ 

t —^0 t t —^0 t 2 t 

where Q{u) = 2 In\^ ~ + ^\^dx - see (22). Using the Taylor series w.r.t t yields 

J'{u)uj= ( W F \ f {u — z)u)dx 

Jn Jq. 

with W = —p,{d — V“u -I- ^). Recall that 

f W • V“a;da: = (-1)” f w‘^div“TVcia;-V(-l)^ f 
Jq Jq . g Jo ’ 


j=o 

n — l 


dx 


n-j-l 




Xi — 1 

c 

312 = 1 

c 

X2—0 


dx-} 


dxi- 


(32) 


(33) 


(34) 


where we note n = 2 for 1 < a < 2. Next consider 2 case studies. 


i). Given u{x)\g^ = bi{x), and 

d(^u(x)+tuj{x)^ 
dn 


du{x) 

dn 


an 


= b 2 {x), since [u{x) + tuj{x))\g^ = (u(a:))|gj.j = bi{x) and 


dQ 


= IQQ = b 2 {x), it suffices to take u G '^g^(D,K). Such a choice ensures ^ 


dQ 

























Jianping Zhang and Ke Chen 


21 



(a) u 



Nvyv' \y\ 

0 o.-t 0.2 0.3 0.4 o.s 

O.s 0.7 O.S O.S 


-- --- 


a_B 0.7 D_S o.s 1 


■n -:- =°lL..lor, 1 .w V. V- 

— —--—=3- 

0 0.1 0.2 0.3 O.-* O.S 

o.s 0.7 O.S O.S 1 


—^ » 



(b) error u — u* 


Figure 5: Comparison I — The slice presentations of four restorations along a diagonal line in Fig. 4(a) 


0,z = 0,1 


0 ^ ^ 




( 33 ) reduces to 


xi=o or 1 


dx^' 


2:2=0 or 1 


= 0, n — j — 1 = 0,1. Hence equation (32) with 


ii). Keep u S IFf (H). Since ^ 


equation (341 can only diminish if 


or 1 ’ dxr^-^ 


x2=o or 1 


7^ 0, the boundary terms in 


[“’*'1 xi=o or 1 

The proof is complete. 


= 0 and 


x2=o or 1 


= 0 ^ •71 = 0,7 = 0,1. 


Remark 5 In imaging applications, the above first set i) of boundary conditions seems not reasonable, because 
one hardly knows a priori what bi, 62 should be. The second set ii) of boundary conditions appears complicated 
which might be simplified as follows. 

From J6'M Section 2 . 3.6 pp. 75 ], ifWi(x) has a sufficient number of continuous derivatives, then 


D^-X’W, 


= 0 for any a G (1, 2) is equivalent to 
xi=o or 1 '7^1 


IFi 


a:i=0 or 1 OXi 


xi=o or 1 

= 0 . 


= 0 (j = 0,1), i.e.. 


xi=o or 1 
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Indeed, if the n-th derivative of u(x) is integrable in [0,1], then Wi 

du{x) 


i{x) 


= 0 and 
xi=o or 1 oxi 


xi=o or 1 


= 0 : 


= 0 is equivalent to 


on the other hand, ^ 


a;i=o or 1 


xi=o or 1 

= 0 (for all k = 0,1,2) are equivalent to ^ 


= 0 and 


xi=o or 1 


dx\^ 


aii=o or 1 


= 0, hence one has 


dxi 


or 1 


= 0. The derivations 0/W2 are similar to those ofWi. 
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